! ----------------------------
!
! Applies 4th order Runge-Kutta to Variables in U
!
! ----------------------------
subroutine rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: K1, K2, K3, K4

  ! Computes the four components
  call deltaU(K1, U             , GAMA, FORCE, A, B, C, m)
  K1 = dt * K1
  call deltaU(K2, U + 0.5d0 * K1, GAMA, FORCE, A, B, C, m)
  K2 = dt * K2
  call deltaU(K3, U + 0.5d0 * K2, GAMA, FORCE, A, B, C, m)
  K3 = dt * K3
  call deltaU(K4, U + K3        , GAMA, FORCE, A, B, C, m)
  K4 = dt * K4

  U = U + (K1 + 2.d0 * K2 + 2.d0 * K3 + K4) / 6.d0

  return
end subroutine rungeKutta4


! ----------------------------
!
! Applies 6th order Runge-Kutta to Variables in U
!
! ----------------------------
subroutine rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: K1, K2, K3, K4, K5, K6, K7

  ! fx_y
  !
  ! f denotes fraction
  ! x - numerator
  ! y - denominator 
  real*8 f1_2
  real*8 f1_3
  real*8 f1_12
  real*8 f1_16
  real*8 f2_3
  real*8 f3_4
  real*8 f3_8
  real*8 f3_16
  real*8 f9_8
  real*8 f9_11
  real*8 f9_44
  real*8 f16_11
  real*8 f18_11
  real*8 f63_44
  
  parameter(f1_2   =  1.d0 /  2.d0)
  parameter(f1_3   =  1.d0 /  3.d0)
  parameter(f1_12  =  1.d0 / 12.d0)
  parameter(f1_16  =  1.d0 / 16.d0)
  parameter(f2_3   =  2.d0 /  3.d0)
  parameter(f3_4   =  3.d0 /  4.d0)
  parameter(f3_8   =  3.d0 /  8.d0)
  parameter(f3_16  =  3.d0 / 16.d0)
  parameter(f9_8   =  9.d0 /  8.d0)
  parameter(f9_11  =  9.d0 / 11.d0)
  parameter(f9_44  =  9.d0 / 44.d0)
  parameter(f16_11 = 16.d0 / 11.d0)
  parameter(f18_11 = 18.d0 / 11.d0)
  parameter(f63_44 = 63.d0 / 44.d0)  
  
  ! ci_j_k_...
  !
  ! c denotes coefficient
  ! i, j, k - index of the corresponding parcels 
  real*8 c1_7
  real*8 c3_4
  real*8 c5_6
  real*8 denominador
  
  parameter(denominador = 120.d0)
  parameter(c1_7 = 11.d0 / denominador)
  parameter(c3_4 = 81.d0 / denominador)
  parameter(c5_6 = 32.d0 / denominador)
  
  
  ! Computes the seven components
  call deltaU(K1, U,                                                           &
              GAMA, FORCE, A, B, C, m)
  K1 = dt * K1
  call deltaU(K2, U + f1_3 * K1,                                               &
              GAMA, FORCE, A, B, C, m)
  K2 = dt * K2
  call deltaU(K3, U + f2_3 * K2,                                               &
              GAMA, FORCE, A, B, C, m)
  K3 = dt * K3
  call deltaU(K4, U + f1_12 * K1 + f1_3 * K2 - f1_12 * K3,                     &
              GAMA, FORCE, A, B, C, m)
  K4 = dt * K4
  call deltaU(K5, U + f9_8 * K2 - f1_16 * K1 - f3_16 * K3 - f3_8 * K4,         &
              GAMA, FORCE, A, B, C, m)
  K5 = dt * K5
  call deltaU(K6, U + f9_8 * K2 - f3_8 * K3 - f3_4 * K4 + f1_2 * K5,           & 
              GAMA, FORCE, A, B, C, m)
  K6 = dt * K6
  call deltaU(K7, U + f9_44*K1 - f9_11*K2 + f63_44*K3 + f18_11*K4 - f16_11*K5, &
              GAMA, FORCE, A, B, C, m)
  K7 = dt * K7
	
  U = U + c1_7 * (K1 + K7) + c3_4 * (K3 + K4) - c5_6 * (K5 + K6)

  return
end subroutine rungeKutta6


! ----------------------------
!
! Initiates Adams-Bashforth Variables To Use With The 4th order Method
!
! ----------------------------
subroutine initAdamBashforth4(U, F0, F1, F2, F3, A, B, C, GAMA, FORCE, dt, m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: F0, F1, F2, F3

  call deltaU(F0, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F1, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F2, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F3, U, GAMA, FORCE, A, B, C, m)

  return
end subroutine initAdamBashforth4


! ----------------------------
!
! Initiates Adams-Bashforth Variables To Use With The 6th order Method
!
! ----------------------------
subroutine initAdamBashforth6(U, F0, F1, F2, F3, F4, F5, A, B, C, GAMA, FORCE, &
                              dt, m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: F0, F1, F2, F3, F4, F5

  call deltaU(F0, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F1, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F2, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F3, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F4, U, GAMA, FORCE, A, B, C, m)
  call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
  call deltaU(F5, U, GAMA, FORCE, A, B, C, m)

  return
end subroutine initAdamBashforth6


! ----------------------------
!
! Applies 4th order Adams-Bashforth to Variables in U
!
! ----------------------------
subroutine AdamBashforth4(U, F0, F1, F2, F3, A, B, C, GAMA, FORCE, dt, m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: F0, F1, F2, F3

  ! ci
  !
  ! c denotes coefficient
  ! i - index of the corresponding parcel
  real*8 c0
  real*8 c1
  real*8 c2
  real*8 c3
  real*8 denominador
  
  parameter(denominador = 24.d0)
  parameter(c0 = 55.d0 / denominador)
  parameter(c1 = 59.d0 / denominador)
  parameter(c2 = 37.d0 / denominador)
  parameter(c3 =  9.d0 / denominador)

  U = U + (c0 * F3 - c1 * F2 + c2 * F1 - c3 * F0) * dt
  call deltaU(F0, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c0 * F0 - c1 * F3 + c2 * F2 - c3 * F1) * dt
  call deltaU(F1, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c0 * F1 - c1 * F0 + c2 * F3 - c3 * F2) * dt
  call deltaU(F2, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c0 * F2 - c1 * F1 + c2 * F0 - c3 * F3) * dt
  call deltaU(F3, U, GAMA, FORCE, A, B, C, m)

  return
end subroutine AdamBashforth4


! ----------------------------
!
! Applies 6th order Adams-Bashforth to Variables in U
!
! ----------------------------
subroutine AdamBashforth6(U, F0, F1, F2, F3, F4, F5, A, B, C, GAMA, FORCE, dt, &
                          m)
  implicit none

  integer*8 m
  real*8 dt
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: F0, F1, F2, F3, F4, F5

  ! ci
  !
  ! c denotes coefficient
  ! i - index of the corresponding parcel
  real*8 c0
  real*8 c1
  real*8 c2
  real*8 c3
  real*8 c4
  real*8 c5
  real*8 denominador
  
  parameter(denominador = 1440.d0)
  parameter(c0 =  475.d0 / denominador)
  parameter(c1 = 2877.d0 / denominador)
  parameter(c2 = 7298.d0 / denominador)
  parameter(c3 = 9982.d0 / denominador)
  parameter(c4 = 7923.d0 / denominador)
  parameter(c5 = 4277.d0 / denominador)

  U = U + (c5 * F5 - c4 * F4 + c3 * F3 - c2 * F2 + c1 * F1 - c0 * F0) * dt
  call deltaU(F0, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c5 * F0 - c4 * F5 + c3 * F4 - c2 * F3 + c1 * F2 - c0 * F1) * dt
  call deltaU(F1, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c5 * F1 - c4 * F0 + c3 * F5 - c2 * F4 + c1 * F3 - c0 * F2) * dt
  call deltaU(F2, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c5 * F2 - c4 * F1 + c3 * F0 - c2 * F5 + c1 * F4 - c0 * F3) * dt
  call deltaU(F3, U, GAMA, FORCE, A, B, C, m)
  
  U = U + (c5 * F3 - c4 * F2 + c3 * F1 - c2 * F0 + c1 * F5 - c0 * F4) * dt
  call deltaU(F4, U, GAMA, FORCE, A, B, C, m)

  U = U + (c5 * F4 - c4 * F3 + c3 * F2 - c2 * F1 + c1 * F0 - c0 * F5) * dt
  call deltaU(F5, U, GAMA, FORCE, A, B, C, m)
    
  return
end subroutine AdamBashforth6

! ----------------------------
!
! Computes the next U according with the previous values
! Returns the result in Udot
!
! ----------------------------
subroutine deltaU(UDOT, U, GAMA, FORCE, A, B, C, m)
  implicit none
  
  integer*8 m
  real*8,     dimension( 1:m  ):: GAMA 
  real*8,     dimension( 1:m  ):: A, B, C
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: UDOT, U

  UDOT(  0) = (0.d0, 0.d0)
  UDOT( -1) = (0.d0, 0.d0)
  UDOT(m+1) = (0.d0, 0.d0)
  UDOT(m+2) = (0.d0, 0.d0)
  UDOT(1:m) = (0.d0, 1.d0) * &
              DCONJG(A * U(2:m+1) * U( 3:m+2) + &
                     B * U(0:m-1) * U( 2:m+1) + &
                     C * U(0:m-1) * U(-1:m-2)) - GAMA * U(1:m) + FORCE

  return
end subroutine deltaU


! ----------------------------
!
! Applies a slaved scheme based in Adams-Bashforth method.
! Does two iterations in time.
!
! ----------------------------
subroutine slavedAdamsBashforth2(U, G0, G1, FORCE, LT, NLT, A, B, C, m)
  implicit none

  integer*8 m
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension( 1:m  ):: FORCE, G0, G1
  real*8,     dimension( 1:m  ):: A, B, C, LT, NLT
  
  U(1:m) = LT * U(1:m) + NLT * (1.5d0 * G1 - 0.5d0 * G0)
  call G(G0, U, FORCE, A, B, C, m)
  U(1:m) = LT * U(1:m) + NLT * (1.5d0 * G0 - 0.5d0 * G1)
  call G(G1, U, FORCE, A, B, C, m)

  return
end subroutine slavedAdamsBashforth2


! ----------------------------
!
! Returns in the variable OUTCOME the result of the function G.
!
! ----------------------------
subroutine G(OUTCOME, U, FORCE, A, B, C, m)
  implicit none

  integer*8 m
  complex*16, dimension(-1: m+2):: U
  complex*16, dimension( 1: m  ):: FORCE, OUTCOME
  real*8,     dimension( 1: m  ):: A, B, C

  OUTCOME = (0.d0, 1.d0) * DCONJG(A * U(2:m+1) * U( 3:m+2) + &
                                  B * U(0:m-1) * U( 2:m+1) + &
                                  C * U(0:m-1) * U(-1:m-2)) + FORCE
  return
end subroutine G


! ----------------------------
!
! Calculates the energy value for the given U.
!
! ----------------------------
real*8 function E(U, m)
  implicit none

  integer*8 n
  integer*8 m
  complex*16, dimension(-1: m+2):: U

  E = 0.d0

  do n = 1, m
    E = E + CDABS(U(n))**2
  end do

  E = E / 2.d0

  return
end function E


! ----------------------------
!
! Calculates the enstrophy value for the given U.
!
! ----------------------------
real*8 function H(U, m)
  implicit none

  integer*8 n
  integer*8 m
  complex*16, dimension(-1: m+2):: U

  H = 0.d0

  do n = 1, m
    H = H + (-2.d0)**n * CDABS(U(n))**2
  end do

  H = H / 2.d0

  return
end function H


! ----------------------------
!
! Computes the structure function for the given velocity matrix.
!
! ----------------------------
subroutine structureFunction(OUTPUT, VELOCITY, order, numberOutputs, m)
  implicit none
  
  integer*8 k
  integer*8 m
  integer*8 n
  integer*4 order
  integer*8 numberOutputs
  real*8    , dimension( 1:m  ):: OUTPUT
  complex*16, dimension( m, numberOutputs):: VELOCITY
     
  OUTPUT = 0.d0

  do k = 1, numberOutputs
    do n = 1, m
      OUTPUT(n) = OUTPUT(n) + CDABS(VELOCITY(n, k))**order
    end do
  end do

  OUTPUT = OUTPUT / numberOutputs
  
  return 
end subroutine structureFunction


! ----------------------------
!
! Computes the third order moment for the given velocity matrix.
!
! ----------------------------
subroutine thirdOrderMoment(OUTPUT, VELOCITY, numberOutputs, m)
  implicit none
  
  integer*8 k
  integer*8 m
  integer*8 n
  integer*8 numberOutputs
  real*8    , dimension(1: m):: OUTPUT
  complex*16, dimension(1: m):: SOMATORY
  complex*16, dimension( m, numberOutputs):: VELOCITY
  
  SOMATORY = (0.d0, 0.d0)
  
  do k = 1, numberOutputs
    do n = 1, m
      SOMATORY(n) = SOMATORY(n) + &
      VELOCITY(n - 1, k) * VELOCITY(n, k) * VELOCITY(n + 1, k) + &
      4.d0 * VELOCITY(n, k) * VELOCITY(n + 1, k) * VELOCITY(n + 2, k)
    end do
  end do
  
  OUTPUT = 0.25d0 * DABS(DIMAG(SOMATORY)) / numberOutputs
  
  return
end subroutine thirdOrderMoment


! ----------------------------
!
! Returns the velocity values of some specific shell.
!
! ----------------------------
subroutine nShellValues(OUTPUT, VELOCITY, n, numberOutputs, m)
  implicit none
  
  integer*8 k
  integer*8 m
  integer*8 n
  integer*8 numberOutputs
  complex*16 VELOCITY(m, numberOutputs)
  complex*16 OUTPUT(numberOutputs)

  do k = 1, numberOutputs
    OUTPUT(k) = VELOCITY(n, k)
  end do
  
  return
end subroutine nShellValues


! ----------------------------
!
! Returns a matrix with the velocity modulus of all shells over all time.
!
! ----------------------------
subroutine shellModulus(OUTPUT, VELOCITY, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 k
  integer*8 n
  integer*8 numberOutputs
  complex*16 VELOCITY(m, numberOutputs)
  real*4 OUTPUT(m, numberOutputs)
  
  do k = 1, numberOutputs
    do n = 1, m
      OUTPUT(n, k) = REAL(CDABS(VELOCITY(n, k)))
    end do
  end do
  
  return
end subroutine shellModulus


! ----------------------------
!
! Reads the values of the matrix Velocity with the values in the fileName file.
!
! ----------------------------
subroutine readVelocity(VELOCITY, fileName, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 k
  integer*8 n
  integer*8 numberOutputs
  real*8 value
  character*(*) fileName
  complex*16 VELOCITY(m, numberOutputs)

  open (10, file = fileName)
  
  do k = 1, numberOutputs
    read (10, *) value
    do n = 1, m
      read (10, *) value
      VELOCITY(n, k) = value
      read (10, *) value
      VELOCITY(n, k) = VELOCITY(n, k) + (0.d0, 1.d0) * value
    end do
  end do  
  close(10)

  return
end subroutine readVelocity
